{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 高斯玻色采样模拟分子振动光谱"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "此研究涉及两个主要领域。其中之一是关于量子信息和量子计算机的研究。我们选择**玻色采样**进行讨论。\n",
    "目前，有多种物理实现方式可以进行玻色子取样，但最常见的是光子解决方案。\n",
    "在该方案中，单光子态被用作输入到干涉仪，在输出端测量光子数。\n",
    "干涉装置相对简单；而制备单光子和测量光子数对仪器的要求较高。\n",
    "因此，如果我们考虑使用相干激光，可以让实验简便一些。\n",
    "由于相干态可以用高斯分布描述，这种类型的玻色子采样被称为高斯玻色采样。\n",
    "\n",
    "\n",
    "另一个主要领域是**光谱学**。我们以分子振动光谱为例，模拟不同带电态之间转变的概率，从而近似确定分子振动谱。\n",
    "这种方法可以使用高斯玻色采样来实现。\n",
    "参考 [Huh](https://doi.org/10.1038/nphoton.2015.153) 的研究工作，我们在基于玻色采样的模拟器上实现该模拟过程。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 高斯玻色采样\n",
    "\n",
    "量子比特是一个抽象的二维量子物理系统。\n",
    "其希尔伯特空间通常可由两个基态 $\\left|0\\right\\rangle$ 和 $\\left|1\\right\\rangle$ 来描述。\n",
    "不同量子比特可以处于叠加态和纠缠状态。\n",
    "我们可将其推广到3维、4维以至有限的 $d$-维空间中。\n",
    "\n",
    "以上内容可以进一步推广到无限维度。\n",
    "玻色子粒子是这种类型的量子物理系统的一个例子，其中在给定状态中可以有无限多个粒子。\n",
    "这样的系统仍然是量子系统，仍具有叠加和纠缠性。\n",
    "在可数无限维度情况下，量子态被称为 qumodes。\n",
    "\n",
    "在真实物理系统中，任何物理量都需要通过多次测量来确定其状态分布。\n",
    "因此量子计算机都是对抽样问题的物理实现。\n",
    "在可数无限维度的情况下，这被称为玻色子采样。\n",
    "\n",
    "在多种玻色采样的物理实现方式中，最简单的仍然是光子方案。\n",
    "量子光学设备在 qumodes 上为幺正算符。常用高斯算符如下：\n",
    "\n",
    "\n",
    "$$\\mathcal{P}_{i}\\left(\\phi\\right)=e^{-i\\phi a_{i}^{\\dagger}a_{i}}$$\n",
    "\n",
    "$$\\mathcal{D}_{i}\\left(\\alpha\\right)=e^{\\alpha a_{i}^{\\dagger}-\\alpha^{*}a_{i}}$$\n",
    "\n",
    "$$\\mathcal{S}_{i}\\left(\\zeta\\right)=e^{\\frac{1}{2}\\left(\\zeta^{*}a_{i}^{2}-\\zeta a_{i}^{\\dagger2}\\right)}$$\n",
    "\n",
    "$$\\mathcal{B}_{ij}\\left(\\theta\\right)=e^{\\theta\\left(a_{i}^{\\dagger}a_{j}-a_{i}a_{j}^{\\dagger}\\right)}$$\n",
    "\n",
    "\n",
    "事实上，一个纯的光子数态，即 $\\left|n\\right\\rangle$, 在技术上很难做到。\n",
    "$\\left|1\\right\\rangle$ 意味着特定量子模式中只有一个光子。\n",
    "精确地产生一个光子是一项巨大挑战。\n",
    "而相干光是每个量子模式的可数无限状态之和，高斯玻色采样以相干光作为输入，是一种更简单的解决方案。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 分子的振动光谱简要描述"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div align=center>\n",
    "\t<img src=\"./fig/Franck_Condon_Diagram.png\" width=\"30%\"/>\t\n",
    "</div>\n",
    "\n",
    "<center>\n",
    "  <a href=\"https://chem.libretexts.org/@api/deki/files/79071/800px-Franck_Condon_Diagram.svg.png?revision=1\">Franck–Condon原理能级图</a>\n",
    "</center>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "光谱学的一个主要问题是检查原子和分子之间电荷状态的跃迁。\n",
    "分子吸收的光频率便是取决于不同电子态之间的允许跃迁。\n",
    "这些电子跃迁可能伴随着分子振动能量的变化。\n",
    "在这种情况下，代表更强烈吸收光的频率的吸收线被称为振动光谱。\n",
    "然而，除了氢原子可以被精确解决外，随着电子数量增加，量子力学方程变得指数级复杂，使得对它们的计算变得不可能。因此必须使用近似方法。\n",
    "\n",
    "如果假设核和电子状态是独立的，对于极其复杂的分子系统，一个很好的近似称为玻恩-奥本海默近似。\n",
    "由于原子核的质量较大，其对电子壳层变化反应更慢，因此在 Franck-Condon 原理中可以将原子核视为恒定。\n",
    "\n",
    "通过用简单的量子谐振子来近似电子的振动，可以将电子的哈密顿算符写成反应坐标的函数。\n",
    "然后，电子势能表面已经可以用抛物线描述。\n",
    "\n",
    "$$\\mathcal{H}=p^{2}+q^{2}$$\n",
    "\n",
    "Duschinsky 提出一种近似方法，不同电子态的简正坐标之间存在线性关系\n",
    "\n",
    "$$q'=U_{D}q+d$$\n",
    "\n",
    "Doktorov 证明，在这些条件下，对于这样的量子谐振子，不同势能面上的量子态 $\\left|\\phi\\right\\rangle$ 和 $\\left|\\phi'\\right\\rangle$ 之间的关系可以用以下方式给出:\n",
    "\n",
    "$$\\left|\\phi'\\right\\rangle =U_{Dokt}\\left|\\phi\\right\\rangle  $$\n",
    "\n",
    "其中，Doktorov 算子定义为\n",
    "\n",
    "$$U_{Dok}=D\\left(d\\right)S^{\\dagger}\\left(\\Omega'\\right)R\\left(U_{Dusch}\\right)S_{\\Omega}\\left(\\Omega\\right)$$\n",
    "\n",
    "- 旋转算符的输入参数是Duschinsky混合矩阵本身\n",
    "\n",
    "- 位移算符D的参数是Duschinsky位移\n",
    "- 压缩算符是从分子的物理特性推导出来的。$\\Omega$ 和 $\\Omega'$ 是分子内原子在发生跃迁前后谐波的角频率\n",
    "\n",
    "上面哈密顿量的本征值问题的解是 Fock 空间中的相干态。高斯算符作用于它们之后，将一个相干态转换为另一个相干态。\n",
    "不同电子态之间的转移概率，即所谓的 Franck-Condon 因子，可以近似表示为\n",
    "\n",
    "$$\n",
    "FCF_{n',n}\t=\\left|\\left\\langle \\phi'|\\phi\\right\\rangle \\right|^{2}\n",
    "\t=\\left|\\left\\langle n'\\right|U_{Dokt}\\left|n\\right\\rangle \\right|^{2}\n",
    "$$\n",
    "\t\n",
    "玻色子采样的计算在于对 qumodes 进行旋转矩阵操作。\n",
    "在我们计算振动光谱时，使用上面的 Doktorov 算符。其过程为\n",
    "\n",
    "1. 制备相干态\n",
    "\n",
    "2. 将相干态传输到一个玻色采样设备，以获得转移概率\n",
    "3. 进行第二次压缩，测量各模式分布，进一步计算FCF\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hessian matrix\n",
    "\n",
    "对于一双原子分子，假设谐振子势能\n",
    "\n",
    "$$V=\\frac{k}{2}(r-r_e)^2$$\n",
    "\n",
    "其中 $r$ 是两个核之间的距离，$r_e$ 是它们的平衡距离。\n",
    "简单起见，可只考虑X轴。我们可以用两个原子的笛卡尔坐标 $x_1$（对应原子1）和 $x_2$（对应原子2）来表示这一点。\n",
    "则，$r = x_2 - x_1$。\n",
    "\n",
    "势能的二阶导数是：\n",
    "\n",
    "$$\\frac{d^{2}V}{dx_{1}^{2}}=k$$\n",
    "\n",
    "$$\\frac{d^{2}V}{dx_{2}^{2}}=k$$\n",
    "\n",
    "$$\\frac{d^{2}V}{dx_{1}dx_{2}}=\\frac{d^{2}V}{dx_{2}dx_{1}}=-k$$\n",
    "\n",
    "\n",
    "Hessian矩阵（仅针对x方向）为：\n",
    "\n",
    "$$\n",
    "H=\\left(\\begin{array}{cc}\n",
    "k & -k\\\\\n",
    "-k & k\n",
    "\\end{array}\\right)\n",
    "$$\n",
    "\n",
    "质量加权的Hessian矩阵是：\n",
    "\n",
    "$$\n",
    "F=\\left(\\begin{array}{cc}\n",
    "\\frac{k}{m_{1}} & -\\frac{k}{\\sqrt{m_{1}m_{2}}}\\\\\n",
    "-\\frac{k}{\\sqrt{m_{1}m_{2}}} & \\frac{k}{m_{2}}\n",
    "\\end{array}\\right)\n",
    "$$\n",
    "\n",
    "如果我们有了系统的Hessian矩阵和质量加权Hessian矩阵，就可以由特征向量定义质量加权简正模式。\n",
    "接着，利用分子的平衡结构坐标数据、质量加权简正坐标等数据可以求解所需高斯算符参数。\n",
    "\n",
    "\n",
    "[Jankowiak](https://pubs.aip.org/aip/jcp/article-abstract/127/23/234101/906357/Vibronic-transitions-in-large-molecular-systems?redirectedFrom=fulltext)的补充材料提供了一些分子的相关数据。\n",
    "我们以甲酸（$1^{1}A'\\rightarrow1^{2}A'$）为例计算分子振动光谱："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import deepquantum as dq\n",
    "import torch\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$1^1A'$ 态及$1^2A'$ 态的平衡结构坐标如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "ri = np.genfromtxt(\"./data/formic_ri.csv\", delimiter=\",\", skip_header=0)[:, np.newaxis]\n",
    "rf = np.genfromtxt(\"./data/formic_rf.csv\", delimiter=\",\", skip_header=0)[:, np.newaxis]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$1^1A'$ 态及 $1^2A'$ 态的质量加权简正坐标如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "li = np.genfromtxt(\"./data/formic_li.csv\", delimiter=\",\", skip_header=0)\n",
    "lf = np.genfromtxt(\"./data/formic_lf.csv\", delimiter=\",\", skip_header=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$1^1A'$ 态及 $1^2A'$ 态简正模式频率如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "omega = np.genfromtxt(\"./data/formic_omega.csv\", delimiter=\",\", skip_header=0)\n",
    "omegap = np.genfromtxt(\"./data/formic_omegap.csv\", delimiter=\",\", skip_header=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "计算中用到的物理学常量"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = 299792458.0  # 光速\n",
    "mu = 1.6605390666 * 10**-27  # 原子质量单位\n",
    "h = 6.62607015 * 10**-34  # 普朗克常数\n",
    "\n",
    "m_c = 12  # 碳原子相对原子质量\n",
    "m_h = 1.007825037  # 氢原子相对原子质量\n",
    "m_o = 15.994914640  # 氧原子相对原子质量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Duschinsky 矩阵 $U$ 及位移矢量 $\\delta$ 计算如下"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.99343181,  0.01440011,  0.01532633,  0.02861045,  0.06378083,\n",
       "         0.07513988, -0.04280796],\n",
       "       [-0.01485231,  0.99314303,  0.07419536,  0.0769291 , -0.03610189,\n",
       "        -0.00248855,  0.01727882],\n",
       "       [-0.01185718, -0.09164895,  0.84227531,  0.17994129, -0.38567948,\n",
       "         0.30738928,  0.08008576],\n",
       "       [ 0.03813492,  0.04087165, -0.34031972, -0.52311845, -0.66785609,\n",
       "         0.38477092,  0.11420873],\n",
       "       [-0.04133251, -0.03419326, -0.40036477,  0.76362651, -0.10356638,\n",
       "         0.48381836,  0.09406589],\n",
       "       [ 0.0907918 , -0.04184572, -0.09066643,  0.31511355, -0.59003401,\n",
       "        -0.719268  ,  0.13037051],\n",
       "       [-0.03245558,  0.00500992, -0.02063661,  0.06935002, -0.20181377,\n",
       "         0.01732396, -0.97588364]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u = []\n",
    "for li_ele in li:\n",
    "    for lf_elf in lf:\n",
    "        u.append(np.sum(li_ele * lf_elf))\n",
    "u = np.array(u[-1::-1]).reshape(7, 7).T\n",
    "u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.22536617],\n",
       "       [ 0.14689208],\n",
       "       [ 1.55989779],\n",
       "       [-0.37838396],\n",
       "       [ 0.45525871],\n",
       "       [-0.34391138],\n",
       "       [ 0.06184607]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "delta = []\n",
    "m = np.diag([m_c, m_c, m_c, m_o, m_o, m_o, m_o, m_o, m_o, m_h, m_h, m_h, m_h, m_h, m_h])\n",
    "for i in range(len(omegap)):\n",
    "    d = lf[i].T @ np.sqrt(m) @ (ri - rf)\n",
    "    l = np.sqrt(h / (4 * np.pi**2 * 100 * omegap[i] * c * mu)) / (10**-10)\n",
    "    delta.append(d / l)\n",
    "delta = np.array(delta[-1::-1])\n",
    "delta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在用于计算振动光谱的GBS算法中，以上这些化学参数足以确定GBS设备的配置。\n",
    "利用这些数据，我们即可计算甲酸分子振动光谱。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "实际上，在可能仅涉及有限数量的光子的情况下需要**非线性相互作用**，第二次挤压操作在光学装置中难以实现。\n",
    "通常，我们需要将两次挤压操作压缩为一次：\n",
    "\n",
    "$$U_{Dok}=R_{C_{L}}S_{\\Sigma}^{\\dagger}R_{C_{R}}D_{2^{-1/2}J^{-1}\\delta}$$\n",
    "\n",
    "其中\n",
    "$$J=\\Omega'U\\Omega^{-1}$$\n",
    "$$\\delta=\\hbar^{-1/2}\\Omega'd$$\n",
    "$$J=C_{L}\\Sigma C_{R}^{t}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "各 GBS 参数计算如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "pre_transition_squeezing = np.sqrt(omega[-1::-1])\n",
    "post_transition_squeezing = np.sqrt(omegap[-1::-1])\n",
    "\n",
    "j_mat = (\n",
    "    np.diag(post_transition_squeezing)\n",
    "    @ u\n",
    "    @ np.linalg.inv(np.diag(pre_transition_squeezing))\n",
    ")\n",
    "\n",
    "cl, lambda_1, cr = np.linalg.svd(j_mat)\n",
    "\n",
    "delta_2 = np.linalg.inv(j_mat) @ delta / np.sqrt(2)\n",
    "delta_2 = delta_2.flatten()\n",
    "lambda_2 = np.log(lambda_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 甲酸分子振动光谱计算如下所示"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "modes = 7  # 简正模式数量\n",
    "cutoff = 3\n",
    "shots = 500000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg baseProfile=\"full\" height=\"6.363636363636363cm\" version=\"1.1\" width=\"12.0cm\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:ev=\"http://www.w3.org/2001/xml-events\" xmlns:xlink=\"http://www.w3.org/1999/xlink\"><defs /><polyline fill=\"none\" points=\"40,30 130,30\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"25\" /><text font-size=\"9\" x=\"83\" y=\"20\">D</text><text font-size=\"7\" x=\"95\" y=\"18\">r =0.056</text><text font-size=\"7\" x=\"95\" y=\"24\">θ =0.0</text><polyline fill=\"none\" points=\"40,60 130,60\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"55\" /><text font-size=\"9\" x=\"83\" y=\"50\">D</text><text font-size=\"7\" x=\"95\" y=\"48\">r =-0.053</text><text font-size=\"7\" x=\"95\" y=\"54\">θ =0.0</text><polyline fill=\"none\" points=\"40,90 130,90\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"85\" /><text font-size=\"9\" x=\"83\" y=\"80\">D</text><text font-size=\"7\" x=\"95\" y=\"78\">r =0.982</text><text font-size=\"7\" x=\"95\" y=\"84\">θ =0.0</text><polyline fill=\"none\" points=\"40,120 130,120\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"115\" /><text font-size=\"9\" x=\"83\" y=\"110\">D</text><text font-size=\"7\" x=\"95\" y=\"108\">r =0.525</text><text font-size=\"7\" x=\"95\" y=\"114\">θ =0.0</text><polyline fill=\"none\" points=\"40,150 130,150\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"145\" /><text font-size=\"9\" x=\"83\" y=\"140\">D</text><text font-size=\"7\" x=\"95\" y=\"138\">r =-0.112</text><text font-size=\"7\" x=\"95\" y=\"144\">θ =0.0</text><polyline fill=\"none\" points=\"40,180 130,180\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"175\" /><text font-size=\"9\" x=\"83\" y=\"170\">D</text><text font-size=\"7\" x=\"95\" y=\"168\">r =0.525</text><text font-size=\"7\" x=\"95\" y=\"174\">θ =0.0</text><polyline fill=\"none\" points=\"40,210 130,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"green\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"82.5\" y=\"205\" /><text font-size=\"9\" x=\"83\" y=\"200\">D</text><text font-size=\"7\" x=\"95\" y=\"198\">r =-0.016</text><text font-size=\"7\" x=\"95\" y=\"204\">θ =0.0</text><polyline fill=\"none\" points=\"130,30 150,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,30 220,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,60 150,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,60 220,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,90 150,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,90 220,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,120 150,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,120 220,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,150 150,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,150 220,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,180 150,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,180 220,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"130,210 150,210\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"200,210 220,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"cadetblue\" height=\"200\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"2\" width=\"50\" x=\"150\" y=\"20\" /><text font-size=\"10\" x=\"171\" y=\"120.0\">U</text><polyline fill=\"none\" points=\"220,30 310,30\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"25\" /><text font-size=\"9\" x=\"263\" y=\"20\">S</text><text font-size=\"7\" x=\"275\" y=\"18\">r =-0.097</text><text font-size=\"7\" x=\"275\" y=\"24\">θ =0.0</text><polyline fill=\"none\" points=\"220,60 310,60\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"55\" /><text font-size=\"9\" x=\"263\" y=\"50\">S</text><text font-size=\"7\" x=\"275\" y=\"48\">r =-0.07</text><text font-size=\"7\" x=\"275\" y=\"54\">θ =0.0</text><polyline fill=\"none\" points=\"220,90 310,90\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"85\" /><text font-size=\"9\" x=\"263\" y=\"80\">S</text><text font-size=\"7\" x=\"275\" y=\"78\">r =-0.021</text><text font-size=\"7\" x=\"275\" y=\"84\">θ =0.0</text><polyline fill=\"none\" points=\"220,120 310,120\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"115\" /><text font-size=\"9\" x=\"263\" y=\"110\">S</text><text font-size=\"7\" x=\"275\" y=\"108\">r =0.06</text><text font-size=\"7\" x=\"275\" y=\"114\">θ =0.0</text><polyline fill=\"none\" points=\"220,150 310,150\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"145\" /><text font-size=\"9\" x=\"263\" y=\"140\">S</text><text font-size=\"7\" x=\"275\" y=\"138\">r =0.075</text><text font-size=\"7\" x=\"275\" y=\"144\">θ =0.0</text><polyline fill=\"none\" points=\"220,180 310,180\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"175\" /><text font-size=\"9\" x=\"263\" y=\"170\">S</text><text font-size=\"7\" x=\"275\" y=\"168\">r =0.112</text><text font-size=\"7\" x=\"275\" y=\"174\">θ =0.0</text><polyline fill=\"none\" points=\"220,210 310,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"royalblue\" height=\"12\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"1.5\" width=\"10\" x=\"262.5\" y=\"205\" /><text font-size=\"9\" x=\"263\" y=\"200\">S</text><text font-size=\"7\" x=\"275\" y=\"198\">r =0.187</text><text font-size=\"7\" x=\"275\" y=\"204\">θ =0.0</text><polyline fill=\"none\" points=\"310,30 330,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,30 400,30\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,60 330,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,60 400,60\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,90 330,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,90 400,90\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,120 330,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,120 400,120\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,150 330,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,150 400,150\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,180 330,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,180 400,180\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"310,210 330,210\" stroke=\"black\" stroke-width=\"2\" /><polyline fill=\"none\" points=\"380,210 400,210\" stroke=\"black\" stroke-width=\"2\" /><rect fill=\"cadetblue\" height=\"200\" rx=\"0\" ry=\"0\" stroke=\"black\" stroke-width=\"2\" width=\"50\" x=\"330\" y=\"20\" /><text font-size=\"10\" x=\"351\" y=\"120.0\">U</text><text font-size=\"12\" x=\"25\" y=\"30\">0</text><text font-size=\"12\" x=\"25\" y=\"60\">1</text><text font-size=\"12\" x=\"25\" y=\"90\">2</text><text font-size=\"12\" x=\"25\" y=\"120\">3</text><text font-size=\"12\" x=\"25\" y=\"150\">4</text><text font-size=\"12\" x=\"25\" y=\"180\">5</text><text font-size=\"12\" x=\"25\" y=\"210\">6</text></svg>"
      ],
      "text/plain": [
       "<svgwrite.drawing.Drawing at 0x309ee5670>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cir = dq.photonic.QumodeCircuit(\n",
    "    nmode=modes,\n",
    "    init_state=\"vac\",\n",
    "    # init_state=init_state,\n",
    "    cutoff=cutoff,\n",
    "    backend=\"gaussian\",\n",
    ")\n",
    "\n",
    "for i in range(modes):\n",
    "    cir.d(wires=[i], r=delta_2[i])\n",
    "\n",
    "cir.any(cr, wires=list(range(modes)))\n",
    "\n",
    "for i in range(modes):\n",
    "    cir.s(wires=[i], r=-lambda_2[i])\n",
    "\n",
    "cir.any(cl, wires=list(range(modes)))\n",
    "\n",
    "state = cir()\n",
    "\n",
    "# 线路可视化\n",
    "cir.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "chain 1: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:05<00:00, 18541.01it/s]\u001b[0m\n",
      "chain 2: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 34781.32it/s]\u001b[0m\n",
      "chain 3: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 34753.39it/s]\u001b[0m\n",
      "chain 4: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 36514.41it/s]\u001b[0m\n",
      "chain 5: 100%|\u001b[32m█████████████████████████\u001b[0m| 99999/99999 [00:02<00:00, 35294.55it/s]\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "sample = cir.measure(shots=shots)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAG4CAYAAAAAMkB2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAA9hAAAPYQGoP6dpAABH7UlEQVR4nO3de1xVdb7/8fdGrmp74yVADJXSvKfmBSmz6ciIZU2UzVEjY4ps6kCplLcpLzU1Nna0tEyOM002p5rK+Y2OqWGEt1JCRVExJZtMLGdDk7K3mgrK9/dHhzXuQEVEWcjr+Xisx8P9/X7WWt/v2hv227X3WjiMMUYAAACwBb+6HgAAAAD+jXAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICN+Nf1ABqS8vJyHThwQFdccYUcDkddDwcAAFSDMUaHDx9WZGSk/Pwu/nktwtkldODAAUVFRdX1MAAAQA3s379fV1111UXfD+HsErriiisk/fjkOp3OOh4NAACoDq/Xq6ioKOt9/GIjnF1CFR9lOp1OwhkAAPXMpfpKEhcEAAAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNEM4AAABsxL+uBwBcLhzPOKpsN9PMJR4JAKA+48wZAACAjRDOAAAAbIRwBgAAYCOEMwAAABshnAEAANgI4QwAAMBGCGcAAAA2QjgDAACwEcIZAACAjRDOAAAAbIRwBgAAYCOEMwAAABshnAEAANhInYazdevW6Y477lBkZKQcDoeWLFlSqWbXrl36xS9+IZfLpSZNmqhv374qLCy0+o8fP66UlBS1aNFCTZs21bBhw1RUVOSzjcLCQg0dOlSNGzdWWFiYxo8fr5MnT/rUrFmzRtdff72CgoLUvn17LVy4sNJY5s2bp3bt2ik4OFgxMTHauHFjrRwHAACACnUazo4ePaoePXpo3rx5Vfb/4x//0IABA9SpUyetWbNG27dv15QpUxQcHGzVjBs3Th988IEWLVqktWvX6sCBA7r77rut/lOnTmno0KEqLS3Vhg0b9Oabb2rhwoWaOnWqVbN3714NHTpUt9xyi/Ly8jR27Fg99NBDWrlypVXz3nvvKS0tTdOmTdOWLVvUo0cPxcfHq7i4+CIcGQAA0FA5jDGmrgchSQ6HQ4sXL1ZCQoLVNmLECAUEBOh///d/q1zH4/Hoyiuv1DvvvKN77rlHkrR792517txZ2dnZ6t+/vz788EPdfvvtOnDggMLDwyVJ6enpmjhxor777jsFBgZq4sSJWr58ufLz8332XVJSooyMDElSTEyM+vbtq1dffVWSVF5erqioKD322GOaNGlStebo9Xrlcrnk8XjkdDrP+xjB3hzPOKpsN9Ns8SMGAKihS/3+bdvvnJWXl2v58uW69tprFR8fr7CwMMXExPh89Jmbm6uysjLFxcVZbZ06dVKbNm2UnZ0tScrOzlb37t2tYCZJ8fHx8nq92rlzp1Vz+jYqaiq2UVpaqtzcXJ8aPz8/xcXFWTVVOXHihLxer88CAABwNrYNZ8XFxTpy5IheeOEFDRkyRB999JHuuusu3X333Vq7dq0kye12KzAwUKGhoT7rhoeHy+12WzWnB7OK/oq+s9V4vV4dO3ZM//rXv3Tq1Kkqayq2UZUZM2bI5XJZS1RU1PkfCAAA0KDYNpyVl5dLku68806NGzdOPXv21KRJk3T77bcrPT29jkdXPZMnT5bH47GW/fv31/WQAACAzdk2nLVs2VL+/v7q0qWLT3vnzp2tqzUjIiJUWlqqkpISn5qioiJFRERYNT+9erPi8blqnE6nQkJC1LJlSzVq1KjKmoptVCUoKEhOp9NnAQAAOBvbhrPAwED17dtXBQUFPu1ffPGF2rZtK0nq3bu3AgIClJWVZfUXFBSosLBQsbGxkqTY2Fjt2LHD56rKzMxMOZ1OK/jFxsb6bKOipmIbgYGB6t27t09NeXm5srKyrBoAAIDa4F+XOz9y5Ii+/PJL6/HevXuVl5en5s2bq02bNho/fryGDx+ugQMH6pZbblFGRoY++OADrVmzRpLkcrmUnJystLQ0NW/eXE6nU4899phiY2PVv39/SdLgwYPVpUsXjRo1SjNnzpTb7dbTTz+tlJQUBQUFSZIeeeQRvfrqq5owYYIefPBBrVq1Su+//76WL19ujS0tLU1JSUnq06eP+vXrp5dffllHjx7VAw88cOkOGAAAuPyZOrR69WojqdKSlJRk1bz++uumffv2Jjg42PTo0cMsWbLEZxvHjh0z//Vf/2WaNWtmGjdubO666y7zz3/+06fm66+/NrfeeqsJCQkxLVu2NE888YQpKyurNJaePXuawMBAc/XVV5s33nij0nhfeeUV06ZNGxMYGGj69etnPvvss/Oar8fjMZKMx+M5r/VQP2i6qlwAAPXbpX7/ts19zhoC7nN2eeM+ZwBweeI+ZwAAAA0Y4QwAAMBGCGcAAAA2QjgDAACwEcIZAACAjRDOAAAAbIRwBgAAYCOEMwAAABshnAEAANgI4QwAAMBGCGcAAAA2QjgDAACwEcIZAACAjRDOAAAAbIRwBgAAYCOEMwAAABshnAEAANgI4QwAAMBGCGcAAAA2QjgDAACwEcIZAACAjRDOAAAAbIRwBgAAYCOEMwAAABshnAEAANgI4QwAAMBGCGcAAAA2QjgDAACwEcIZAACAjRDOAAAAbKROw9m6det0xx13KDIyUg6HQ0uWLDlj7SOPPCKHw6GXX37Zp/3gwYNKTEyU0+lUaGiokpOTdeTIEZ+a7du366abblJwcLCioqI0c+bMSttftGiROnXqpODgYHXv3l0rVqzw6TfGaOrUqWrVqpVCQkIUFxenPXv21HjuAAAAVanTcHb06FH16NFD8+bNO2vd4sWL9dlnnykyMrJSX2Jionbu3KnMzEwtW7ZM69at08MPP2z1e71eDR48WG3btlVubq5efPFFTZ8+XQsWLLBqNmzYoJEjRyo5OVlbt25VQkKCEhISlJ+fb9XMnDlTc+fOVXp6unJyctSkSRPFx8fr+PHjtXAkAAAA/o+xCUlm8eLFldq/+eYb07p1a5Ofn2/atm1rXnrpJavv888/N5LMpk2brLYPP/zQOBwO8+233xpjjHnttddMs2bNzIkTJ6yaiRMnmo4dO1qP//M//9MMHTrUZ78xMTHm17/+tTHGmPLychMREWFefPFFq7+kpMQEBQWZv/zlL9Weo8fjMZKMx+Op9jqoPzRdVS4AgPrtUr9/2/o7Z+Xl5Ro1apTGjx+vrl27VurPzs5WaGio+vTpY7XFxcXJz89POTk5Vs3AgQMVGBho1cTHx6ugoECHDh2yauLi4ny2HR8fr+zsbEnS3r175Xa7fWpcLpdiYmKsmqqcOHFCXq/XZwEAADgbW4ez3//+9/L399fjjz9eZb/b7VZYWJhPm7+/v5o3by63223VhIeH+9RUPD5Xzen9p69XVU1VZsyYIZfLZS1RUVFnnS8AAIBtw1lubq7mzJmjhQsXyuFw1PVwamTy5MnyeDzWsn///roeEgAAsDnbhrNPPvlExcXFatOmjfz9/eXv7699+/bpiSeeULt27SRJERERKi4u9lnv5MmTOnjwoCIiIqyaoqIin5qKx+eqOb3/9PWqqqlKUFCQnE6nzwIAAHA2tg1no0aN0vbt25WXl2ctkZGRGj9+vFauXClJio2NVUlJiXJzc631Vq1apfLycsXExFg169atU1lZmVWTmZmpjh07qlmzZlZNVlaWz/4zMzMVGxsrSYqOjlZERIRPjdfrVU5OjlUDAABQG/zrcudHjhzRl19+aT3eu3ev8vLy1Lx5c7Vp00YtWrTwqQ8ICFBERIQ6duwoSercubOGDBmi0aNHKz09XWVlZUpNTdWIESOs227ce++9euaZZ5ScnKyJEycqPz9fc+bM0UsvvWRtd8yYMbr55ps1a9YsDR06VO+++642b95s3W7D4XBo7Nixeu6559ShQwdFR0drypQpioyMVEJCwkU+SgAAoCGp03C2efNm3XLLLdbjtLQ0SVJSUpIWLlxYrW28/fbbSk1N1aBBg+Tn56dhw4Zp7ty5Vr/L5dJHH32klJQU9e7dWy1bttTUqVN97oV2ww036J133tHTTz+t3/zmN+rQoYOWLFmibt26WTUTJkzQ0aNH9fDDD6ukpEQDBgxQRkaGgoODL/AoAAAA/JvDGGPqehANhdfrlcvlksfj4ftnlyHHM1VfuGKm8SMGAPXZpX7/tu13zgAAABoiwhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEbqNJytW7dOd9xxhyIjI+VwOLRkyRKrr6ysTBMnTlT37t3VpEkTRUZG6v7779eBAwd8tnHw4EElJibK6XQqNDRUycnJOnLkiE/N9u3bddNNNyk4OFhRUVGaOXNmpbEsWrRInTp1UnBwsLp3764VK1b49BtjNHXqVLVq1UohISGKi4vTnj17au9gAAAAqI7D2dGjR9WjRw/NmzevUt8PP/ygLVu2aMqUKdqyZYv+9re/qaCgQL/4xS986hITE7Vz505lZmZq2bJlWrdunR5++GGr3+v1avDgwWrbtq1yc3P14osvavr06VqwYIFVs2HDBo0cOVLJycnaunWrEhISlJCQoPz8fKtm5syZmjt3rtLT05WTk6MmTZooPj5ex48fvwhHBgAANFQOY4yp60FIksPh0OLFi5WQkHDGmk2bNqlfv37at2+f2rRpo127dqlLly7atGmT+vTpI0nKyMjQbbfdpm+++UaRkZGaP3++nnrqKbndbgUGBkqSJk2apCVLlmj37t2SpOHDh+vo0aNatmyZta/+/furZ8+eSk9PlzFGkZGReuKJJ/Tkk09Kkjwej8LDw7Vw4UKNGDGiWnP0er1yuVzyeDxyOp01OUywMcczjirbzTRb/IgBAGroUr9/16vvnHk8HjkcDoWGhkqSsrOzFRoaagUzSYqLi5Ofn59ycnKsmoEDB1rBTJLi4+NVUFCgQ4cOWTVxcXE++4qPj1d2drYkae/evXK73T41LpdLMTExVk1VTpw4Ia/X67MAAACcTb0JZ8ePH9fEiRM1cuRIK7W63W6FhYX51Pn7+6t58+Zyu91WTXh4uE9NxeNz1Zzef/p6VdVUZcaMGXK5XNYSFRV1XnMGAAANT70IZ2VlZfrP//xPGWM0f/78uh5OtU2ePFkej8da9u/fX9dDAgAANudf1wM4l4pgtm/fPq1atcrns96IiAgVFxf71J88eVIHDx5URESEVVNUVORTU/H4XDWn91e0tWrVyqemZ8+eZxx7UFCQgoKCzme6AACggbP1mbOKYLZnzx59/PHHatGihU9/bGysSkpKlJuba7WtWrVK5eXliomJsWrWrVunsrIyqyYzM1MdO3ZUs2bNrJqsrCyfbWdmZio2NlaSFB0drYiICJ8ar9ernJwcqwYAAKA21Gk4O3LkiPLy8pSXlyfpxy/e5+XlqbCwUGVlZbrnnnu0efNmvf322zp16pTcbrfcbrdKS0slSZ07d9aQIUM0evRobdy4UevXr1dqaqpGjBihyMhISdK9996rwMBAJScna+fOnXrvvfc0Z84cpaWlWeMYM2aMMjIyNGvWLO3evVvTp0/X5s2blZqaKunHK0nHjh2r5557TkuXLtWOHTt0//33KzIy8qxXlwIAAJw3U4dWr15tJFVakpKSzN69e6vsk2RWr15tbeP77783I0eONE2bNjVOp9M88MAD5vDhwz772bZtmxkwYIAJCgoyrVu3Ni+88EKlsbz//vvm2muvNYGBgaZr165m+fLlPv3l5eVmypQpJjw83AQFBZlBgwaZgoKC85qvx+MxkozH4zmv9VA/aLqqXAAA9dulfv+2zX3OGgLuc3Z54z5nAHB54j5nAAAADRjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEYIZwAAADZCOAMAALAR/7oeAOo3/tg3AAC1izNnAAAANkI4AwAAsBHCGQAAgI0QzgAAAGyEcAYAAGAjhDMAAAAbIZwBAADYCOEMAADARghnAAAANkI4AwAAsBHCGQAAgI0QzgAAAGyEcAYAAGAjhDMAAAAbIZwBAADYCOEMAADARghnAAAANkI4AwAAsJE6DWfr1q3THXfcocjISDkcDi1ZssSn3xijqVOnqlWrVgoJCVFcXJz27NnjU3Pw4EElJibK6XQqNDRUycnJOnLkiE/N9u3bddNNNyk4OFhRUVGaOXNmpbEsWrRInTp1UnBwsLp3764VK1ac91gAAAAuVJ2Gs6NHj6pHjx6aN29elf0zZ87U3LlzlZ6erpycHDVp0kTx8fE6fvy4VZOYmKidO3cqMzNTy5Yt07p16/Twww9b/V6vV4MHD1bbtm2Vm5urF198UdOnT9eCBQusmg0bNmjkyJFKTk7W1q1blZCQoISEBOXn55/XWAAAAC6Uwxhj6noQkuRwOLR48WIlJCRI+vFMVWRkpJ544gk9+eSTkiSPx6Pw8HAtXLhQI0aM0K5du9SlSxdt2rRJffr0kSRlZGTotttu0zfffKPIyEjNnz9fTz31lNxutwIDAyVJkyZN0pIlS7R7925J0vDhw3X06FEtW7bMGk///v3Vs2dPpaenV2ss1eH1euVyueTxeOR0OmvluNU1xzOOKtvNNFu8rC4pjgUAXJ4u9fu3bb9ztnfvXrndbsXFxVltLpdLMTExys7OliRlZ2crNDTUCmaSFBcXJz8/P+Xk5Fg1AwcOtIKZJMXHx6ugoECHDh2yak7fT0VNxX6qM5aqnDhxQl6v12cBAAA4G9uGM7fbLUkKDw/3aQ8PD7f63G63wsLCfPr9/f3VvHlzn5qqtnH6Ps5Uc3r/ucZSlRkzZsjlcllLVFTUOWYNAAAaOtuGs8vB5MmT5fF4rGX//v11PSQAAGBztg1nERERkqSioiKf9qKiIqsvIiJCxcXFPv0nT57UwYMHfWqq2sbp+zhTzen95xpLVYKCguR0On0WAACAs6lRONuyZYt27NhhPf773/+uhIQE/eY3v1FpaWmtDCw6OloRERHKysqy2rxer3JychQbGytJio2NVUlJiXJzc62aVatWqby8XDExMVbNunXrVFZWZtVkZmaqY8eOatasmVVz+n4qair2U52xAAAA1IYahbNf//rX+uKLLyRJX331lUaMGKHGjRtr0aJFmjBhQrW3c+TIEeXl5SkvL0/Sj1+8z8vLU2FhoRwOh8aOHavnnntOS5cu1Y4dO3T//fcrMjLSuqKzc+fOGjJkiEaPHq2NGzdq/fr1Sk1N1YgRIxQZGSlJuvfeexUYGKjk5GTt3LlT7733nubMmaO0tDRrHGPGjFFGRoZmzZql3bt3a/r06dq8ebNSU1MlqVpjAQAAqBWmBpxOp/nyyy+NMca88MILZvDgwcYYYz799FNz1VVXVXs7q1evNpIqLUlJScYYY8rLy82UKVNMeHi4CQoKMoMGDTIFBQU+2/j+++/NyJEjTdOmTY3T6TQPPPCAOXz4sE/Ntm3bzIABA0xQUJBp3bq1eeGFFyqN5f333zfXXnutCQwMNF27djXLly/36a/OWM7F4/EYScbj8ZzXenam6apyaYg4FgBwebrU7981us+Z0+lUbm6uOnTooJ///Oe6/fbbNWbMGBUWFqpjx446duxYrQbIywX3Obu8cSwA4PJUL+5z1qdPHz333HP63//9X61du1ZDhw6V9OPHkj+93QQAAACqr0bh7KWXXtKWLVuUmpqqp556Su3bt5ck/fWvf9UNN9xQqwMEAABoSPxrslKPHj18rtas8OKLL8rfv0abBAAAgGp45uzqq6/W999/X6n9+PHjuvbaay94UAAAAA1VjcLZ119/rVOnTlVqP3HihL755psLHhQAAEBDdV6fQS5dutT698qVK+VyuazHp06dUlZWlqKjo2tvdAAAAA3MeYWzihuuOhwOJSUl+fQFBASoXbt2mjVrVq0NDgAAoKE5r3BWXl4u6cc/Z7Rp0ya1bNnyogwKAACgoarRpZV79+6t7XEAAABANQxnkpSVlaWsrCwVFxdbZ9Qq/OlPf7rggQEAADRENQpnzzzzjJ599ln16dNHrVq1ksNR9Z+tAQAAwPmpUThLT0/XwoULNWrUqNoeDwAAQINWo/uclZaW8meaAAAALoIahbOHHnpI77zzTm2PBQAAoMGr0ceax48f14IFC/Txxx/ruuuuU0BAgE//7Nmza2VwAAAADU2Nwtn27dvVs2dPSVJ+fr5PHxcHAAAA1FyNwtnq1atrexwAAABQDb9zBgAAgIujRmfObrnllrN+fLlq1aoaDwgAAKAhq1E4q/i+WYWysjLl5eUpPz+/0h9EBwAAQPXVKJy99NJLVbZPnz5dR44cuaABAQAANGS1+p2z++67j7+rCQAAcAFqNZxlZ2crODi4NjcJAADQoNToY827777b57ExRv/85z+1efNmTZkypVYGBgAA0BDVKJy5XC6fx35+furYsaOeffZZDR48uFYGBgAA0BDVKJy98cYbtT0OAAAAqIbhrEJubq527dolSeratat69epVK4MCAABoqGoUzoqLizVixAitWbNGoaGhkqSSkhLdcsstevfdd3XllVfW5hgBAAAajBpdrfnYY4/p8OHD2rlzpw4ePKiDBw8qPz9fXq9Xjz/+eG2PEQAAoMGo0ZmzjIwMffzxx+rcubPV1qVLF82bN48LAgAAAC5Ajc6clZeXKyAgoFJ7QECAysvLL3hQAAAADVWNwtl//Md/aMyYMTpw4IDV9u2332rcuHEaNGhQrQ3u1KlTmjJliqKjoxUSEqJrrrlGv/3tb2WMsWqMMZo6dapatWqlkJAQxcXFac+ePT7bOXjwoBITE+V0OhUaGqrk5ORKf2Zq+/btuummmxQcHKyoqCjNnDmz0ngWLVqkTp06KTg4WN27d9eKFStqba4AAABSDcPZq6++Kq/Xq3bt2umaa67RNddco+joaHm9Xr3yyiu1Nrjf//73mj9/vl599VXt2rVLv//97zVz5kyffcycOVNz585Venq6cnJy1KRJE8XHx+v48eNWTWJionbu3KnMzEwtW7ZM69at08MPP2z1e71eDR48WG3btlVubq5efPFFTZ8+XQsWLLBqNmzYoJEjRyo5OVlbt25VQkKCEhISlJ+fX2vzBQAAcJjTT0OdB2OMPv74Y+3evVuS1LlzZ8XFxdXq4G6//XaFh4fr9ddft9qGDRumkJAQvfXWWzLGKDIyUk888YSefPJJSZLH41F4eLgWLlyoESNGaNeuXerSpYs2bdqkPn36SPrxO3O33XabvvnmG0VGRmr+/Pl66qmn5Ha7FRgYKEmaNGmSlixZYs1v+PDhOnr0qJYtW2aNpX///urZs6fS09OrNR+v1yuXyyWPxyOn01krx6iuOZ5xVNluptXoZVWvcSwA4PJ0qd+/z+vM2apVq9SlSxd5vV45HA79/Oc/12OPPabHHntMffv2VdeuXfXJJ5/U2uBuuOEGZWVl6YsvvpAkbdu2TZ9++qluvfVWSdLevXvldrt9QqHL5VJMTIyys7Ml/fj3PkNDQ61gJklxcXHy8/NTTk6OVTNw4EArmElSfHy8CgoKdOjQIavmp+EzPj7e2k9VTpw4Ia/X67MAAACczXmFs5dfflmjR4+uMjW6XC79+te/1uzZs2ttcJMmTdKIESPUqVMnBQQEqFevXho7dqwSExMlSW63W5IUHh7us154eLjV53a7FRYW5tPv7++v5s2b+9RUtY3T93Gmmor+qsyYMUMul8taoqKizmv+AACg4TmvcLZt2zYNGTLkjP2DBw9Wbm7uBQ+qwvvvv6+3335b77zzjrZs2aI333xT//3f/60333yz1vZxMU2ePFkej8da9u/fX9dDAgAANnde9zkrKiqq8hYa1sb8/fXdd99d8KAqjB8/3jp7Jkndu3fXvn37NGPGDCUlJSkiIsIaV6tWrXzG2bNnT0lSRESEiouLfbZ78uRJHTx40Fo/IiJCRUVFPjUVj89VU9FflaCgIAUFBZ3vtAEAQAN2XmfOWrdufdarE7dv3+4Tki7UDz/8ID8/3yE2atTIupdadHS0IiIilJWVZfV7vV7l5OQoNjZWkhQbG6uSkhKfM3qrVq1SeXm5YmJirJp169aprKzMqsnMzFTHjh3VrFkzq+b0/VTUVOwHAACgNpxXOLvttts0ZcoUn9tUVDh27JimTZum22+/vdYGd8cdd+j555/X8uXL9fXXX2vx4sWaPXu27rrrLkmSw+HQ2LFj9dxzz2np0qXasWOH7r//fkVGRiohIUHSj1eRDhkyRKNHj9bGjRu1fv16paamasSIEYqMjJQk3XvvvQoMDFRycrJ27typ9957T3PmzFFaWpo1ljFjxigjI0OzZs3S7t27NX36dG3evFmpqam1Nl8AAIDzupVGUVGRrr/+ejVq1Eipqanq2LGjJGn37t2aN2+eTp06pS1btlT64nxNHT58WFOmTNHixYtVXFysyMhIjRw5UlOnTrWurDTGaNq0aVqwYIFKSko0YMAAvfbaa7r22mut7Rw8eFCpqan64IMP5Ofnp2HDhmnu3Llq2rSpVbN9+3alpKRo06ZNatmypR577DFNnDjRZzyLFi3S008/ra+//lodOnTQzJkzddttt1V7PtxK4/LGsQCAy9Olfv8+7/uc7du3T48++qhWrlxp3anf4XAoPj5e8+bNU3R09EUZ6OWAcHZ541gAwOXpUr9/n/cfPm/btq1WrFihQ4cO6csvv5QxRh06dLC+mwUAAICaO+9wVqFZs2bq27dvbY4FAACgwavR39YEAADAxUE4AwAAsBHCGQAAgI0QzgAAAGyEcAYAAGAjhDMAAAAbIZwBAADYSI3vcwYAZ8JfSwCAmuPMGQAAgI0QzgAAAGyEcAYAAGAjhDMAAAAbIZwBAADYCOEMAADARghnAAAANkI4AwAAsBHCGQAAgI0QzgAAAGyEcAYAAGAjhDMAAAAbIZwBAADYCOEMAADARghnAAAANkI4AwAAsBHCGQAAgI0QzgAAAGyEcAYAAGAjhDMAAAAbIZwBAADYiO3D2bfffqv77rtPLVq0UEhIiLp3767Nmzdb/cYYTZ06Va1atVJISIji4uK0Z88en20cPHhQiYmJcjqdCg0NVXJyso4cOeJTs337dt10000KDg5WVFSUZs6cWWksixYtUqdOnRQcHKzu3btrxYoVF2fSAACgwbJ1ODt06JBuvPFGBQQE6MMPP9Tnn3+uWbNmqVmzZlbNzJkzNXfuXKWnpysnJ0dNmjRRfHy8jh8/btUkJiZq586dyszM1LJly7Ru3To9/PDDVr/X69XgwYPVtm1b5ebm6sUXX9T06dO1YMECq2bDhg0aOXKkkpOTtXXrViUkJCghIUH5+fmX5mAAAIAGwWGMMXU9iDOZNGmS1q9fr08++aTKfmOMIiMj9cQTT+jJJ5+UJHk8HoWHh2vhwoUaMWKEdu3apS5dumjTpk3q06ePJCkjI0O33XabvvnmG0VGRmr+/Pl66qmn5Ha7FRgYaO17yZIl2r17tyRp+PDhOnr0qJYtW2btv3///urZs6fS09OrNR+v1yuXyyWPxyOn01nj42InjmccVbababZ9WV00HIt/41gAuJxc6vdvW585W7p0qfr06aNf/vKXCgsLU69evfSHP/zB6t+7d6/cbrfi4uKsNpfLpZiYGGVnZ0uSsrOzFRoaagUzSYqLi5Ofn59ycnKsmoEDB1rBTJLi4+NVUFCgQ4cOWTWn76eipmI/VTlx4oS8Xq/PAgAAcDa2DmdfffWV5s+frw4dOmjlypV69NFH9fjjj+vNN9+UJLndbklSeHi4z3rh4eFWn9vtVlhYmE+/v7+/mjdv7lNT1TZO38eZair6qzJjxgy5XC5riYqKOq/5AwCAhsfW4ay8vFzXX3+9fve736lXr156+OGHNXr06Gp/jFjXJk+eLI/HYy379++v6yEBAACbs3U4a9Wqlbp06eLT1rlzZxUWFkqSIiIiJElFRUU+NUVFRVZfRESEiouLffpPnjypgwcP+tRUtY3T93Gmmor+qgQFBcnpdPosAAAAZ2PrcHbjjTeqoKDAp+2LL75Q27ZtJUnR0dGKiIhQVlaW1e/1epWTk6PY2FhJUmxsrEpKSpSbm2vVrFq1SuXl5YqJibFq1q1bp7KyMqsmMzNTHTt2tK4MjY2N9dlPRU3FfgAAAGqDrcPZuHHj9Nlnn+l3v/udvvzyS73zzjtasGCBUlJSJEkOh0Njx47Vc889p6VLl2rHjh26//77FRkZqYSEBEk/nmkbMmSIRo8erY0bN2r9+vVKTU3ViBEjFBkZKUm69957FRgYqOTkZO3cuVPvvfee5syZo7S0NGssY8aMUUZGhmbNmqXdu3dr+vTp2rx5s1JTUy/5cQEAAJcv/7oewNn07dtXixcv1uTJk/Xss88qOjpaL7/8shITE62aCRMm6OjRo3r44YdVUlKiAQMGKCMjQ8HBwVbN22+/rdTUVA0aNEh+fn4aNmyY5s6da/W7XC599NFHSklJUe/evdWyZUtNnTrV515oN9xwg9555x09/fTT+s1vfqMOHTpoyZIl6tat26U5GAAAoEGw9X3OLjfc5+zyxrH4N44FgMsJ9zkDAABowAhnAAAANkI4AwAAsBHCGQAAgI0QzgAAAGyEcAYAAGAjhDMAAAAbIZwBAADYCOEMAADARghnAAAANkI4AwAAsBHCGQAAgI0QzgAAAGyEcAYAAGAjhDMAAAAbIZwBAADYCOEMAADARghnAAAANkI4AwAAsBHCGQAAgI0QzgAAAGyEcAYAAGAjhDMAAAAbIZwBAADYCOEMAADARghnAAAANkI4AwAAsBHCGQAAgI0QzgAAAGyEcAYAAGAj9SqcvfDCC3I4HBo7dqzVdvz4caWkpKhFixZq2rSphg0bpqKiIp/1CgsLNXToUDVu3FhhYWEaP368Tp486VOzZs0aXX/99QoKClL79u21cOHCSvufN2+e2rVrp+DgYMXExGjjxo0XY5oAAKABqzfhbNOmTfqf//kfXXfddT7t48aN0wcffKBFixZp7dq1OnDggO6++26r/9SpUxo6dKhKS0u1YcMGvfnmm1q4cKGmTp1q1ezdu1dDhw7VLbfcory8PI0dO1YPPfSQVq5cadW89957SktL07Rp07Rlyxb16NFD8fHxKi4uvviTBwAADUa9CGdHjhxRYmKi/vCHP6hZs2ZWu8fj0euvv67Zs2frP/7jP9S7d2+98cYb2rBhgz777DNJ0kcffaTPP/9cb731lnr27Klbb71Vv/3tbzVv3jyVlpZKktLT0xUdHa1Zs2apc+fOSk1N1T333KOXXnrJ2tfs2bM1evRoPfDAA+rSpYvS09PVuHFj/elPf7q0BwMAAFzW6kU4S0lJ0dChQxUXF+fTnpubq7KyMp/2Tp06qU2bNsrOzpYkZWdnq3v37goPD7dq4uPj5fV6tXPnTqvmp9uOj4+3tlFaWqrc3FyfGj8/P8XFxVk1VTlx4oS8Xq/PAgAAcDb+dT2Ac3n33Xe1ZcsWbdq0qVKf2+1WYGCgQkNDfdrDw8PldrutmtODWUV/Rd/Zarxer44dO6ZDhw7p1KlTVdbs3r37jGOfMWOGnnnmmepNFAAAQDY/c7Z//36NGTNGb7/9toKDg+t6OOdt8uTJ8ng81rJ///66HhIAALA5W4ez3NxcFRcX6/rrr5e/v7/8/f21du1azZ07V/7+/goPD1dpaalKSkp81isqKlJERIQkKSIiotLVmxWPz1XjdDoVEhKili1bqlGjRlXWVGyjKkFBQXI6nT4LAADA2dg6nA0aNEg7duxQXl6etfTp00eJiYnWvwMCApSVlWWtU1BQoMLCQsXGxkqSYmNjtWPHDp+rKjMzM+V0OtWlSxer5vRtVNRUbCMwMFC9e/f2qSkvL1dWVpZVAwAAUBts/Z2zK664Qt26dfNpa9KkiVq0aGG1JycnKy0tTc2bN5fT6dRjjz2m2NhY9e/fX5I0ePBgdenSRaNGjdLMmTPldrv19NNPKyUlRUFBQZKkRx55RK+++qomTJigBx98UKtWrdL777+v5cuXW/tNS0tTUlKS+vTpo379+unll1/W0aNH9cADD1yiowEAABoCW4ez6njppZfk5+enYcOG6cSJE4qPj9drr71m9Tdq1EjLli3To48+qtjYWDVp0kRJSUl69tlnrZro6GgtX75c48aN05w5c3TVVVfpj3/8o+Lj462a4cOH67vvvtPUqVPldrvVs2dPZWRkVLpIAAAA4EI4jDGmrgfRUHi9XrlcLnk8nsvm+2eOZxxVtptpDe9lxbH4N44FgMvJpX7/tvV3zgAAABoawhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEZsH85mzJihvn376oorrlBYWJgSEhJUUFDgU3P8+HGlpKSoRYsWatq0qYYNG6aioiKfmsLCQg0dOlSNGzdWWFiYxo8fr5MnT/rUrFmzRtdff72CgoLUvn17LVy4sNJ45s2bp3bt2ik4OFgxMTHauHFjrc8ZAAA0XLYPZ2vXrlVKSoo+++wzZWZmqqysTIMHD9bRo0etmnHjxumDDz7QokWLtHbtWh04cEB333231X/q1CkNHTpUpaWl2rBhg958800tXLhQU6dOtWr27t2roUOH6pZbblFeXp7Gjh2rhx56SCtXrrRq3nvvPaWlpWnatGnasmWLevToofj4eBUXF1+agwEAAC57DmOMqetBnI/vvvtOYWFhWrt2rQYOHCiPx6Mrr7xS77zzju655x5J0u7du9W5c2dlZ2erf//++vDDD3X77bfrwIEDCg8PlySlp6dr4sSJ+u677xQYGKiJEydq+fLlys/Pt/Y1YsQIlZSUKCMjQ5IUExOjvn376tVXX5UklZeXKyoqSo899pgmTZp0zrF7vV65XC55PB45nc7aPjR1wvGMo8p2M61evaxqBcfi3zgWAC4nl/r92/Znzn7K4/FIkpo3by5Jys3NVVlZmeLi4qyaTp06qU2bNsrOzpYkZWdnq3v37lYwk6T4+Hh5vV7t3LnTqjl9GxU1FdsoLS1Vbm6uT42fn5/i4uKsmp86ceKEvF6vzwJIP4aXny4AAEj1LJyVl5dr7NixuvHGG9WtWzdJktvtVmBgoEJDQ31qw8PD5Xa7rZrTg1lFf0Xf2Wq8Xq+OHTumf/3rXzp16lSVNRXb+KkZM2bI5XJZS1RUVM0mDgAAGox6Fc5SUlKUn5+vd999t66HUi2TJ0+Wx+Oxlv3799f1kAAAgM351/UAqis1NVXLli3TunXrdNVVV1ntERERKi0tVUlJic/Zs6KiIkVERFg1P72qsuJqztNrfnqFZ1FRkZxOp0JCQtSoUSM1atSoypqKbfxUUFCQgoKCajZhAADQINn+zJkxRqmpqVq8eLFWrVql6Ohon/7evXsrICBAWVlZVltBQYEKCwsVGxsrSYqNjdWOHTt8rqrMzMyU0+lUly5drJrTt1FRU7GNwMBA9e7d26emvLxcWVlZVg0AAMCFsv2Zs5SUFL3zzjv6+9//riuuuML6fpfL5VJISIhcLpeSk5OVlpam5s2by+l06rHHHlNsbKz69+8vSRo8eLC6dOmiUaNGaebMmXK73Xr66aeVkpJindl65JFH9Oqrr2rChAl68MEHtWrVKr3//vtavny5NZa0tDQlJSWpT58+6tevn15++WUdPXpUDzzwwKU/MAAA4LJk+3A2f/58SdLPfvYzn/Y33nhDv/rVryRJL730kvz8/DRs2DCdOHFC8fHxeu2116zaRo0aadmyZXr00UcVGxurJk2aKCkpSc8++6xVEx0dreXLl2vcuHGaM2eOrrrqKv3xj39UfHy8VTN8+HB99913mjp1qtxut3r27KmMjIxKFwkAAADUVL27z1l9xn3OLm/ncyyqqr2cjhmvCwCXE+5zBgAA0IARzgAAAGzE9t85A4DL/WNgADgdZ84AAABshDNnAOoMFw4AQGWcOQMAALARwhkAAICNEM4AAABshHAGAABgI1wQAAD1CLcVAS5/nDkDAACwEc6cAUADwa1LgPqBM2cAAAA2QjgDAACwEcIZAACAjfCdM8CmuCoPABomzpwBAADYCOEMAADARghnAAAANkI4AwAAsBEuCABwWePGqwDqG86cAQAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNcLUmANgAV5VeWhxv2BlnzgAAAGyEM2dAPcUfRgeAyxNnzgAAAGyEM2cAAKDWcFb/whHOztO8efP04osvyu12q0ePHnrllVfUr1+/uh4W0KDZ6c3ATmPBpcPzjtpEODsP7733ntLS0pSenq6YmBi9/PLLio+PV0FBgcLCwup6eA0CvwDt70xXwTVkpx+TS/V6vRTPw7l+Hvl5xdlwxeyZEc7Ow+zZszV69Gg98MADkqT09HQtX75cf/rTnzRp0qQ6Hh3qEoHk0nM847gov8Tr4xvGxToWlws7/HwSVO3PTj/7hLNqKi0tVW5uriZPnmy1+fn5KS4uTtnZ2VWuc+LECZ04ccJ67PF4JEler/fiDvZSOl5180WbYxX7s8XxPMNxkM4wvurM41w11TwWrhmuSm2eyZ5z1lRVVy3ncyyq+/o5V10V/Y7JZ39DruivNMcLfU2fZf4V+63yuJ5jv2d6jn5aV91xnHXd/3PO10UNXqM/fV6qOhbVec2elws4FtX+2ailn9efOtfzfsbxVFOt/uxXqO5cz/C8VPWz65nsqf3XRTXH4/V6rfEbc4nOfJtLtad67sCBA2rdurU2bNig2NhYq33ChAlau3atcnJyKq0zffp0PfPMM5dymAAA4CL5xz/+oauvvvqi74czZxfR5MmTlZaWZj0uKSlR27ZtVVhYKJerev8bqi+8Xq+ioqK0f/9+OZ3Ouh5OrWJu9dPlPDfp8p4fc6ufLue5eTwetWnTRs2bN78k+yOcVVPLli3VqFEjFRUV+bQXFRUpIiKiynWCgoIUFBRUqd3lcl12L9wKTqeTudVDzK3+upznx9zqp8t5bn5+l+b2sNyEtpoCAwPVu3dvZWVlWW3l5eXKysry+ZgTAADgQnDm7DykpaUpKSlJffr0Ub9+/fTyyy/r6NGj1tWbAAAAF4pwdh6GDx+u7777TlOnTpXb7VbPnj2VkZGh8PDwaq0fFBSkadOmVflRZ33H3Oon5lZ/Xc7zY271E3OrPVytCQAAYCN85wwAAMBGCGcAAAA2QjgDAACwEcIZAACAjRDOasHzzz+vG264QY0bN1ZoaGiVNYWFhRo6dKgaN26ssLAwjR8/XidPnvSpWbNmja6//noFBQWpffv2WrhwYaXtzJs3T+3atVNwcLBiYmK0cePGizCj82fXcVVYt26d7rjjDkVGRsrhcGjJkiU+/cYYTZ06Va1atVJISIji4uK0Z88en5qDBw8qMTFRTqdToaGhSk5O1pEjR3xqtm/frptuuknBwcGKiorSzJkzL/bUNGPGDPXt21dXXHGFwsLClJCQoIKCAp+a48ePKyUlRS1atFDTpk01bNiwSjdUrq3XaG2aP3++rrvuOuumlrGxsfrwww/r/byq8sILL8jhcGjs2LFWW32d3/Tp0+VwOHyWTp061ft5Vfj222913333qUWLFgoJCVH37t21efNmq78+/z5p165dpefO4XAoJSVFUv197k6dOqUpU6YoOjpaISEhuuaaa/Tb3/7W529l2up5M7hgU6dONbNnzzZpaWnG5XJV6j958qTp1q2biYuLM1u3bjUrVqwwLVu2NJMnT7ZqvvrqK9O4cWOTlpZmPv/8c/PKK6+YRo0amYyMDKvm3XffNYGBgeZPf/qT2blzpxk9erQJDQ01RUVFl2KaZ2TXcZ1uxYoV5qmnnjJ/+9vfjCSzePFin/4XXnjBuFwus2TJErNt2zbzi1/8wkRHR5tjx45ZNUOGDDE9evQwn332mfnkk09M+/btzciRI61+j8djwsPDTWJiosnPzzd/+ctfTEhIiPmf//mfizq3+Ph488Ybb5j8/HyTl5dnbrvtNtOmTRtz5MgRq+aRRx4xUVFRJisry2zevNn079/f3HDDDVZ/bb1Ga9vSpUvN8uXLzRdffGEKCgrMb37zGxMQEGDy8/Pr9bx+auPGjaZdu3bmuuuuM2PGjLHa6+v8pk2bZrp27Wr++c9/Wst3331X7+dljDEHDx40bdu2Nb/61a9MTk6O+eqrr8zKlSvNl19+adXU598nxcXFPs9bZmamkWRWr15tjKm/z93zzz9vWrRoYZYtW2b27t1rFi1aZJo2bWrmzJlj1djpeSOc1aI33nijynC2YsUK4+fnZ9xut9U2f/5843Q6zYkTJ4wxxkyYMMF07drVZ73hw4eb+Ph463G/fv1MSkqK9fjUqVMmMjLSzJgxo5Zncn7sOq4z+Wk4Ky8vNxEREebFF1+02kpKSkxQUJD5y1/+Yowx5vPPPzeSzKZNm6yaDz/80DgcDvPtt98aY4x57bXXTLNmzazn1BhjJk6caDp27HiRZ+SruLjYSDJr1641xvw4l4CAALNo0SKrZteuXUaSyc7ONsbU3mv0UmjWrJn54x//eNnM6/Dhw6ZDhw4mMzPT3HzzzVY4q8/zmzZtmunRo0eVffV5Xsb8+DM9YMCAM/Zfbr9PxowZY6655hpTXl5er5+7oUOHmgcffNCn7e677zaJiYnGGPs9b3yseQlkZ2ere/fuPjerjY+Pl9fr1c6dO62auLg4n/Xi4+OVnZ0tSSotLVVubq5PjZ+fn+Li4qyaumDXcZ2PvXv3yu12+8zB5XIpJibGmkN2drZCQ0PVp08fqyYuLk5+fn7KycmxagYOHKjAwECrJj4+XgUFBTp06NAlms2Pf6BXkvUHenNzc1VWVuYzv06dOqlNmzY+87vQ1+jFdurUKb377rs6evSoYmNjL5t5paSkaOjQoZXGUN/nt2fPHkVGRurqq69WYmKiCgsLL4t5LV26VH369NEvf/lLhYWFqVevXvrDH/5g9V9Ov09KS0v11ltv6cEHH5TD4ajXz90NN9ygrKwsffHFF5Kkbdu26dNPP9Wtt94qyX7PG+HsEnC73ZX+ikDFY7fbfdYar9erY8eO6V//+pdOnTpVZU3FNuqCXcd1PirGebY5uN1uhYWF+fT7+/urefPm53wOT9/HxVZeXq6xY8fqxhtvVLdu3ax9BwYGVvo+5E/nd6Gv0Ytlx44datq0qYKCgvTII49o8eLF6tKlS72flyS9++672rJli2bMmFGprz7PLyYmRgsXLlRGRobmz5+vvXv36qabbtLhw4fr9bwk6auvvtL8+fPVoUMHrVy5Uo8++qgef/xxvfnmmz7juxx+nyxZskQlJSX61a9+Ze23vj53kyZN0ogRI9SpUycFBASoV69eGjt2rBITE33GZpfnjT/fdAaTJk3S73//+7PW7Nq1y+dLrkBdS0lJUX5+vj799NO6Hkqt6dixo/Ly8uTxePTXv/5VSUlJWrt2bV0P64Lt379fY8aMUWZmpoKDg+t6OLWq4myEJF133XWKiYlR27Zt9f777yskJKQOR3bhysvL1adPH/3ud7+TJPXq1Uv5+flKT09XUlJSHY+udr3++uu69dZbFRkZWddDuWDvv/++3n77bb3zzjvq2rWr8vLyNHbsWEVGRtryeePM2Rk88cQT2rVr11mXq6++ulrbioiIqHQ1S8XjiIiIs9Y4nU6FhISoZcuWatSoUZU1FduoC3Yd1/moGOfZ5hAREaHi4mKf/pMnT+rgwYPnfA5P38fFlJqaqmXLlmn16tW66qqrrPaIiAiVlpaqpKSk0tjOZ+zneo1eLIGBgWrfvr169+6tGTNmqEePHpozZ069n1dubq6Ki4t1/fXXy9/fX/7+/lq7dq3mzp0rf39/hYeH1+v5nS40NFTXXnutvvzyy3r/vLVq1UpdunTxaevcubP1se3l8vtk3759+vjjj/XQQw9ZbfX5uRs/frx19qx79+4aNWqUxo0bZ521ttvzRjg7gyuvvFKdOnU663L6Z8pnExsbqx07dvg8qZmZmXI6ndYPeWxsrLKysnzWy8zMVGxsrKQf36B69+7tU1NeXq6srCyrpi7YdVznIzo6WhERET5z8Hq9ysnJseYQGxurkpIS5ebmWjWrVq1SeXm5YmJirJp169aprKzMqsnMzFTHjh3VrFmzizZ+Y4xSU1O1ePFirVq1StHR0T79vXv3VkBAgM/8CgoKVFhY6DO/C32NXirl5eU6ceJEvZ/XoEGDtGPHDuXl5VlLnz59lJiYaP27Ps/vdEeOHNE//vEPtWrVqt4/bzfeeGOlW9V88cUXatu2raT6//ukwhtvvKGwsDANHTrUaqvPz90PP/wgPz/fyNOoUSOVl5dLsuHzdl6XD6BK+/btM1u3bjXPPPOMadq0qdm6davZunWrOXz4sDHm35cWDx482OTl5ZmMjAxz5ZVXVnlp8fjx482uXbvMvHnzqryVRlBQkFm4cKH5/PPPzcMPP2xCQ0N9roqpC3Yd1+kOHz5sPS+SzOzZs83WrVvNvn37jDE/XkIdGhpq/v73v5vt27ebO++8s8pLqHv16mVycnLMp59+ajp06OBzCXVJSYkJDw83o0aNMvn5+ebdd981jRs3vuiXvj/66KPG5XKZNWvW+FwC/8MPP1g1jzzyiGnTpo1ZtWqV2bx5s4mNjTWxsbFWf229RmvbpEmTzNq1a83evXvN9u3bzaRJk4zD4TAfffRRvZ7XmZx+taYx9Xd+TzzxhFmzZo3Zu3evWb9+vYmLizMtW7Y0xcXF9Xpexvx42xN/f3/z/PPPmz179pi3337bNG7c2Lz11ltWTX3+fWLMj1fct2nTxkycOLFSX3197pKSkkzr1q2tW2n87W9/My1btjQTJkywauz0vBHOakFSUpKRVGmpuC+MMcZ8/fXX5tZbbzUhISGmZcuW5oknnjBlZWU+21m9erXp2bOnCQwMNFdffbV54403Ku3rlVdeMW3atDGBgYGmX79+5rPPPrvIs6seu46rwurVq6t8jpKSkowxP15GPWXKFBMeHm6CgoLMoEGDTEFBgc82vv/+ezNy5EjTtGlT43Q6zQMPPGAF8Arbtm0zAwYMMEFBQaZ169bmhRdeuOhzq2peknxeP8eOHTP/9V//ZZo1a2YaN25s7rrrLvPPf/7TZzu19RqtTQ8++KBp27atCQwMNFdeeaUZNGiQFczq87zO5KfhrL7Ob/jw4aZVq1YmMDDQtG7d2gwfPtznPmD1dV4VPvjgA9OtWzcTFBRkOnXqZBYsWODTX59/nxhjzMqVK42kSmM2pv4+d16v14wZM8a0adPGBAcHm6uvvto89dRTPre8sNPz5jDmtNvjAgAAoE7xnTMAAAAbIZwBAADYCOEMAADARghnAAAANkI4AwAAsBHCGQAAgI0QzgAAAGyEcAYAAGAjhDMAAAAbIZwBQANw1113qVmzZrrnnnvqeigAzoFwBgANwJgxY/TnP/+5rocBoBoIZwBwEX3//fcKCwvT119/Xafj+NnPfqYrrriiyr4RI0Zo1qxZl3hEAM6EcAbAFn71q1/J4XBUWoYMGVLXQ7sgzz//vO688061a9eurodyRk8//bSef/55eTyeuh4KAEn+dT0AAKgwZMgQvfHGGz5tQUFBF3WfpaWlCgwMvCjb/uGHH/T6669r5cqVF2X7p+vZs6dOnjxZqf2jjz5SZGTkWdft1q2brrnmGr311ltKSUm5WEMEUE2cOQNgG0FBQYqIiPBZmjVrJunHj+Uef/xxTZgwQc2bN1dERISmT5/us355eblmzJih6OhohYSEqEePHvrrX//qU/Ozn/1MqampGjt2rFq2bKn4+HgdPnxYiYmJatKkiVq1aqWXXnpJP/vZzzR27FhJ0p///Ge1aNFCJ06c8NlWQkKCRo0adcb5rFixQkFBQerfv79Pe2FhoZKSkhQeHm6N89NPP5Ukff3113I4HPp//+//aeDAgQoJCVHfvn1VWFioTz75RP3791fjxo01aNAglZSUWNvMy8tTfn5+peVcwazCHXfcoXfffbdatQAuLsIZgHrjzTffVJMmTZSTk6OZM2fq2WefVWZmptU/Y8YM/fnPf1Z6erp27typcePG6b777tPatWsrbScwMFDr169Xenq60tLStH79ei1dulSZmZn65JNPtGXLFqv+l7/8pU6dOqWlS5dabcXFxVq+fLkefPDBM473k08+Ue/evX3a9u3bp379+unYsWNaunSptm/frtTUVDmdTknStm3bJEnz58/X7373O23YsEFFRUW677779MILL+jVV1/V6tWrtW3btkpnGS9Ev379tHHjxkoBFEAdMABgA0lJSaZRo0amSZMmPsvzzz9vjDHm5ptvNgMGDPBZp2/fvmbixInGGGOOHz9uGjdubDZs2OBTk5ycbEaOHGk9vvnmm02vXr2sx16v1wQEBJhFixZZbSUlJaZx48ZmzJgxVtujjz5qbr31VuvxrFmzzNVXX23Ky8vPOKc777zTPPjggz5tt956q7nzzjvPuM706dNN8+bNzb/+9S+r7b777jPt2rUzR48etdqGDBliJkyYcMbt/NSgQYNMy5YtTUhIiGndunWl47Rt2zYjyXz99dfV3iaAi4PvnAGwjVtuuUXz58/3aWvevLn17+uuu86nr1WrViouLpYkffnll/rhhx/085//3KemtLRUvXr18mk7/WzWV199pbKyMvXr189qc7lc6tixo886o0ePVt++ffXtt9+qdevWWrhwoXURw5kcO3ZMwcHB1uN9+/bpww8/1NatW8+4zrZt23TXXXepRYsWVlthYaGGDx+uxo0b+7TdeeedZ9zOT3388cdn7Q8JCZH04/fkANQtwhkA22jSpInat29/xv6AgACfxw6HQ+Xl5ZKkI0eOSJKWL1+u1q1b+9T99KKCJk2anPfYevXqpR49eujPf/6zBg8erJ07d2r58uVnXadly5Y6dOiQ9TgvL0+BgYHq2bPnGdfJy8vT5MmTfdq2bdumcePGWY+PHz+ugoIC9ejR47zncSYHDx6UJF155ZW1tk0ANUM4A3BZ6NKli4KCglRYWKibb7652utdffXVCggI0KZNm9SmTRtJksfj0RdffKGBAwf61D700EN6+eWX9e233youLk5RUVFn3XavXr301ltvWY8DAgJ08uRJ/fDDDz5nwSp4vV59/fXXPmf69u7dK4/H49O2Y8cOGWPUvXv3as/zXPLz83XVVVepZcuWtbZNADVDOANgGydOnJDb7fZp8/f3r1ZguOKKK/Tkk09q3LhxKi8v14ABA+TxeLR+/Xo5nU4lJSWdcb2kpCSNHz9ezZs3V1hYmKZNmyY/P79KH1nee++9evLJJ/WHP/yhWnfbj4+P1+TJk3Xo0CE1a9ZMMTExcrlcevTRRzVp0iQZY7Ru3ToNGjRIHTp00LZt29SoUSN169bN2kZeXp6aN2+utm3b+rRdc801atq06TnHUF2ffPKJBg8eXGvbA1BzXK0JwDYyMjLUqlUrn2XAgAHVXv+3v/2tpkyZohkzZqhz584aMmSIli9frujo6LOuN3v2bMXGxur2229XXFycbrzxRnXu3Nnn+2LSj99FGzZsmJo2baqEhIRzjqd79+66/vrr9f7770uSWrRooQ8++EB79uxR3759NWDAAC1dulRhYWGSfvz4smPHjj773bZtW6XvzG3btq1WP9I8fvy4lixZotGjR9faNgHUnMMYY+p6EABgJ0ePHlXr1q01a9YsJScn+/QNGjRIXbt21dy5c6u1reXLl2v8+PHKz8+Xn589/z88f/58LV68WB999FFdDwWA+FgTALR161bt3r1b/fr1k8fj0bPPPitJPldDHjp0SGvWrNGaNWv02muvVXvbQ4cO1Z49e/Ttt9+e8ztqdSUgIECvvPJKXQ8DwP/hzBmABm/r1q166KGHVFBQoMDAQPXu3VuzZ8/2+cJ9u3btdOjQIU2ZMkVPPvlkHY4WwOWOcAYAAGAj9vwCBAAAQANFOAMAALARwhkAAICNEM4AAABshHAGAABgI4QzAAAAGyGcAQAA2AjhDAAAwEYIZwAAADZCOAMAALARwhkAAICNEM4AAABshHAGAABgI/8fiSmwZIi0x3AAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wave_number = []\n",
    "counts = []\n",
    "for ele in sample.items():\n",
    "    # print(ele[0].state)\n",
    "    wave_number.append(torch.sum(ele[0].state * omegap))\n",
    "    counts.append(ele[1])\n",
    "\n",
    "plt.bar(wave_number, counts, width=100, color=\"g\")\n",
    "\n",
    "plt.xlabel(r\"Energy $(cm^{-1})$\")\n",
    "plt.ylabel(r\"Counts\")\n",
    "plt.xlim(-1000, 8000)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 参考结果\n",
    "\n",
    "\n",
    "<div align=center>\n",
    "\t<img src=\"./fig/formic_spec.png\" width=\"60%\"/>\t\n",
    "</div>\n",
    "\n",
    "<center>\n",
    "  <a href=\"https://strawberryfields.ai/photonics/apps/run_tutorial_vibronic.html\">甲酸分子振动光谱图</a>\n",
    "</center>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "piquasso",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
